Shallow-curved shell structure with geometric nonlinearities

Finite Element Model
Generate model
Dynamical system setup
Linear Modal analysis and SSM setup
Forced response surface using SSMs

Finite Element Model

FEM used in the following reference:

Jain, S., & Tiso, P. (2018). Simulation-free hyper-reduction for geometrically nonlinear structural dynamics: a quadratic manifold lifting approach. Journal of Computational and Nonlinear Dynamics, 13(7), 071003. <https://doi.org/10.1115/1.4040021>

Finite element code taken from the following package:

Jain, S., Marconi, J., Tiso P. (2020). YetAnotherFEcode (Version v1.1). Zenodo. <http://doi.org/10.5281/zenodo.4011282>

Discretization parameters

clc; clear all;
nDiscretization = 10; % Discretization parameter (#DOFs is proportional to the square of this number)
epsilon = 10; % converge at order 5

Generate model

[M,C,K,fnl,f_0,outdof] = build_model(nDiscretization);
n = length(M); % number of degrees of freedom
disp(['Number of degrees of freedom = ' num2str(n)])
disp(['Phase space dimensionality = ' num2str(2*n)])

Building FE model
Assembling M,C,K matrices
Applying boundary conditions
Solving undamped eigenvalue problem

ans = 

  Patch (Deformed Mesh) with properties:

    FaceColor: 'interp'
    FaceAlpha: 1
    EdgeColor: [0 0 0]
    LineStyle: '-'
        Faces: [400×3 double]
     Vertices: [1200×3 double]

  Use GET to show all properties

Using Rayleigh damping
Assembling external force vector
Getting nonlinearity coefficients
Loaded tensors from storage
Total time spent on model assembly = 00:00:13
Number of degrees of freedom = 1320
Phase space dimensionality = 2640

Dynamical system setup

We consider the forced system

$\mathbf{M}\ddot{\mathbf{x}}+\mathbf{C}\dot{\mathbf{x}}+\mathbf{K}\mathbf{x}+\mathbf{f}(\mathbf{x},\dot{\mathbf{x}})=\epsilon\mathbf{f}^{ext}(\mathbf{\Omega}t),$

which can be written in the first-order form as

$\mathbf{B}\dot{\mathbf{z}}=\mathbf{Az}+\mathbf{F}(\mathbf{z})+\epsilon\mathbf{F}^{ext}(\mathbf{\phi}),\\\dot{\mathbf{\phi}} =\mathbf{\Omega}$

where

$\mathbf{z}=\left[\begin{array}{c}\mathbf{x}\\\dot{\mathbf{x}}\end{array}\right],\quad\mathbf{A}=\left[\begin{array}{cc}-\mathbf{K} & \mathbf{0}\\\mathbf{0} & \mathbf{M}\end{array}\right],\mathbf{B}=\left[\begin{array}{cc}\mathbf{C} & \mathbf{M}\\\mathbf{M} & \mathbf{0}\end{array}\right],\quad\quad\mathbf{F}(\mathbf{z})=\left[\begin{array}{c}\mathbf{-\mathbf{f}(\mathbf{x},\dot{\mathbf{x}})}\\\mathbf{0}\end{array}\right],\quad\mathbf{F}^{ext}(\mathbf{z},\mathbf{\phi})=\left[\begin{array}{c}\mathbf{f}^{ext}(\mathbf{\phi})\\\mathbf{0}\end{array}\right]$ .

DSorder = 2;
DS = DynamicalSystem(DSorder);
set(DS,'M',M,'C',C,'K',K,'fnl',fnl);
set(DS.Options,'Emax',5,'Nmax',10,'notation','multiindex')
% set(DS.Options,'Emax',5,'Nmax',10,'notation','tensor')

We assume periodic forcing of the form

$\mathbf{f}^{ext}(\phi) = \mathbf{f}_0\cos(\phi)=\frac{\mathbf{f}_0}{2}e^{i\phi} + \frac{\mathbf{f}_0}{2}e^{-i\phi}$

Fourier coefficients of Forcing

kappas = [-1; 1];
coeffs = 1e-2*[f_0 f_0]/2;
DS.add_forcing(coeffs, kappas,epsilon);

Linear Modal analysis and SSM setup

[V,D,W] = DS.linear_spectral_analysis();

Due to high-dimensionality, we compute only the first 5 eigenvalues with the smallest magnitude. These would also be used to compute the spectral quotients
Assuming a proportional damping hypthesis with symmetric matrices
modal damping ratio for 1 mode is 2.000000e-03
modal damping ratio for 2 mode is 2.000000e-03
modal damping ratio for 3 mode is 2.102524e-03
modal damping ratio for 4 mode is 2.141338e-03
modal damping ratio for 5 mode is 2.369557e-03
the left eigenvectors may be incorrect in case of asymmetry of matrices

 The first 10 nonzero eigenvalues are given as 
   1.0e+02 *

  -0.0030 + 1.4922i
  -0.0030 - 1.4922i
  -0.0060 + 2.9878i
  -0.0060 - 2.9878i
  -0.0071 + 3.3973i
  -0.0071 - 3.3973i
  -0.0076 + 3.5356i
  -0.0076 - 3.5356i
  -0.0101 + 4.2617i
  -0.0101 - 4.2617i

Choose Master subspace (perform resonance analysis)

S = SSM(DS);
set(S.Options, 'reltol', 0.1,'notation','multiindex')
% set(S.Options, 'reltol', 0.1,'notation','tensor')
masterModes = [1,2,3,4];
mFreq = [1 2];
S.choose_E(masterModes);

No (near) outer resonances detected in the (truncated) spectrum
sigma_out = 3
(near) inner resonance detected for the following combination of master eigenvalues
     0     1     1     0
     1     0     1     1
     2     1     0     0
     1     0     0     1
     0     1     1     1
     1     2     0     0
     2     0     0     0
     0     0     2     1
     1     1     1     0
     0     2     0     0
     0     0     1     2
     1     1     0     1

These are in resonance with the follwing eigenvalues of the master subspace
   1.0e+02 *

  -0.0030 + 1.4922i
  -0.0030 + 1.4922i
  -0.0030 + 1.4922i
  -0.0030 - 1.4922i
  -0.0030 - 1.4922i
  -0.0030 - 1.4922i
  -0.0060 + 2.9878i
  -0.0060 + 2.9878i
  -0.0060 + 2.9878i
  -0.0060 - 2.9878i
  -0.0060 - 2.9878i
  -0.0060 - 2.9878i

sigma_in = 3

Forced response surface using SSMs

Obtaining forced response curve in reduced-polar coordinate

set(S.Options, 'reltol', 0.5,'IRtol',0.05,'notation', 'multiindex')

Choose computational domain

set(S.FRCOptions,'coordinates','cartesian','initialSolver','fsolve');
set(S.contOptions, 'h_min', 1e-2,'h_max',2,'PtMX',30);
set(S.Options, 'contribNonAuto',false, 'COMPtype', 'first');
order  = 5;
omega0 = imag(S.E.spectrum(1));
omegaRange = omega0*[0.92 1.07];
epsRange   = [0.01 1]*epsilon;
parRange   = {omegaRange,epsRange};
% we take a previous stored solution as the initial point for 2D
% continuation
set(S.contOptions,'atlasdim',2, 'R', 2,'R_max',100,'R_min',1e-4,'almax',15,'PtMX',10000);
sol = ep_read_solution('','isol-3.ep',3);
z0 = [sol.x(1)*cos(sol.x(2)); sol.x(1)*cos(sol.x(2)); sol.x(3)*cos(sol.x(4)); sol.x(3)*cos(sol.x(4))];
p0 = [sol.p(1); sol.p(2)*100];
% we rescale observables such that all of them have comparable magnitude
optdof = outdof;
scale_states = [1000 1000 1000 1000];
scale_obs    = 1e4*ones(2*numel(outdof)+1,1);
set(S.FRSOptions, 'calFRS', true);
startfrs = tic;
S.extract_FRS('shell_cart',masterModes,order,mFreq,parRange,outdof,optdof,scale_states,scale_obs,{p0,z0});
timings.FRS = toc(startfrs);

(near) outer resonance detected for the following combination of master eigenvalues
     2     0     0     0
     0     0     2     1
     1     1     1     0
     0     2     0     0
     0     0     1     2
     1     1     0     1
     2     0     0     0
     0     0     2     1
     1     1     1     0
     0     2     0     0
     0     0     1     2
     1     1     0     1
     1     0     1     0
     0     1     2     0
     3     0     0     0
     0     1     0     1
     0     3     0     0
     1     0     0     2

These are in resonance with the follwing eigenvalues of the slave subspace
   1.0e+02 *

  -0.0071 + 3.3973i
  -0.0071 + 3.3973i
  -0.0071 + 3.3973i
  -0.0071 - 3.3973i
  -0.0071 - 3.3973i
  -0.0071 - 3.3973i
  -0.0076 + 3.5356i
  -0.0076 + 3.5356i
  -0.0076 + 3.5356i
  -0.0076 - 3.5356i
  -0.0076 - 3.5356i
  -0.0076 - 3.5356i
  -0.0101 + 4.2617i
  -0.0101 + 4.2617i
  -0.0101 + 4.2617i
  -0.0101 - 4.2617i
  -0.0101 - 4.2617i
  -0.0101 - 4.2617i

sigma_out = 3
(near) inner resonance detected for the following combination of master eigenvalues
     0     1     1     0
     1     0     1     1
     2     1     0     0
     1     0     0     1
     0     1     1     1
     1     2     0     0
     2     0     0     0
     0     0     2     1
     1     1     1     0
     0     2     0     0
     0     0     1     2
     1     1     0     1

These are in resonance with the follwing eigenvalues of the master subspace
   1.0e+02 *

  -0.0030 + 1.4922i
  -0.0030 + 1.4922i
  -0.0030 + 1.4922i
  -0.0030 - 1.4922i
  -0.0030 - 1.4922i
  -0.0030 - 1.4922i
  -0.0060 + 2.9878i
  -0.0060 + 2.9878i
  -0.0060 + 2.9878i
  -0.0060 - 2.9878i
  -0.0060 - 2.9878i
  -0.0060 - 2.9878i

sigma_in = 3
Due to (near) outer resonance, the exisitence of the manifold is questionable and the underlying computation may suffer.
Attempting manifold computation
Manifold computation time at order 2 = 00:00:00
Estimated memory usage at order  2 = 4.41E+00 MB
Manifold computation time at order 3 = 00:00:01
Estimated memory usage at order  3 = 8.92E+00 MB
Manifold computation time at order 4 = 00:00:04
Estimated memory usage at order  4 = 1.88E+01 MB
Manifold computation time at order 5 = 00:00:10
Estimated memory usage at order  5 = 3.48E+01 MB

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.


 Run='shell_cart.FRSep': Continue equilibria along FRS.

    STEP   DAMPING               NORMS              COMPUTATION TIMES
  IT SIT     GAMMA     ||d||     ||f||     ||U||   F(x)  DF(x)  SOLVE
   0                          4.52e-17  2.02e+02    0.0    0.0    0.0

 STEP      TIME        ||U||  LABEL  TYPE            om          eps         Rez1         Rez2         Imz1         Imz2          amp       l2z990       l2z660     linfz990     linfz660         rho1         rho2
    0  00:00:00   2.0175e+02      1  EP      1.4226e+02   1.0000e+01  -3.1474e+00  -9.7058e-02   1.3715e-01   1.2571e-02   1.1721e+00   1.1710e+00   5.0963e-02   1.6669e+00   7.7923e-02   3.1504e+00   9.7869e-02
    2  00:00:02   2.0384e+02      2  EP      1.4372e+02   1.0000e+01  -4.0185e+00  -1.9772e-01   2.2624e-01   3.2868e-02   1.4997e+00   1.4963e+00   1.0129e-01   2.1329e+00   1.5301e-01   4.0248e+00   2.0043e-01
   10  00:00:06   2.0089e+02      3          1.4177e+02   8.5260e+00  -2.4975e+00  -5.7214e-02   1.0090e-01   6.8888e-03   9.2947e-01   9.2897e-01   3.0313e-02   1.3208e+00   4.6523e-02   2.4995e+00   5.7627e-02
   20  00:00:11   2.0694e+02      4          1.4599e+02   7.6347e+00  -5.6366e+00  -6.2986e-01   6.2923e-01   2.0098e-01   2.1340e+00   2.1101e+00   3.1859e-01   3.0116e+00   4.7089e-01   5.6716e+00   6.6115e-01
   26  00:00:14   1.9867e+02      5  EP      1.4010e+02   1.0000e+01  -2.3901e+00  -4.3063e-02   7.8637e-02   4.2261e-03   8.8907e-01   8.8876e-01   2.3545e-02   1.2636e+00   3.6571e-02   2.3914e+00   4.3270e-02
   30  00:00:17   2.0600e+02      6  EP      1.4517e+02   1.0000e+01  -5.7215e+00  -5.3169e-01   4.8021e-01   1.2856e-01   2.1526e+00   2.1358e+00   2.6819e-01   3.0518e+00   3.9980e-01   5.7416e+00   5.4702e-01
   36  00:00:21   1.9457e+02      7  EP      1.3729e+02   8.8005e+00  -1.6059e+00  -1.4956e-02   4.0228e-02   1.1198e-03   5.9704e-01   5.9698e-01   8.6276e-03   8.4742e-01   1.3652e-02   1.6064e+00   1.4998e-02
   40  00:00:23   2.0815e+02      8          1.4702e+02   4.2507e+00  -4.7506e+00  -6.0704e-01   8.4551e-01   3.1322e-01   1.8238e+00   1.7952e+00   3.2203e-01   2.5543e+00   4.7055e-01   4.8252e+00   6.8309e-01
  
   . 
   . 
   . 
   
 6980  01:21:42   2.2472e+02   1148  EP      1.5397e+02   1.0000e+01  -3.4021e+00  -1.3799e+01   2.9689e+01  -4.0970e+00   1.2461e+01   1.0976e+01   5.8997e+00   1.7099e+01   9.0790e+00   2.9884e+01   1.4394e+01
 6981  01:21:44   2.2499e+02   1149  EP      1.5407e+02   1.0000e+01   2.6852e-01  -1.4460e+01   3.0253e+01  -5.9974e-01   1.2592e+01   1.1113e+01   5.9213e+00   1.7318e+01   9.1221e+00   3.0254e+01   1.4472e+01
 6982  01:21:47   2.2502e+02   1150  EP      1.5408e+02   1.0000e+01   9.3710e-01  -1.4464e+01   3.0264e+01   4.0937e-02   1.2598e+01   1.1122e+01   5.9166e+00   1.7339e+01   9.1158e+00   3.0278e+01   1.4464e+01
 6984  01:21:50   2.2490e+02   1151  EP      1.5403e+02   1.0000e+01  -1.4408e+00  -1.4289e+01   3.0098e+01  -2.2374e+00   1.2552e+01   1.1068e+01   5.9214e+00   1.7247e+01   9.1182e+00   3.0133e+01   1.4463e+01
 6986  01:21:53   2.2495e+02   1152  EP      1.5405e+02   1.0000e+01  -6.0715e-01  -1.4401e+01   3.0196e+01  -1.4394e+00   1.2576e+01   1.1094e+01   5.9235e+00   1.7293e+01   9.1238e+00   3.0202e+01   1.4473e+01

Plot results

runid = 'shell_cart.FRSep';
% FRS in physical coordinates
titles = {'$||w_\mathrm{A}||_{\mathcal{L}^2}$'
    '$||w_\mathrm{B}||_{\mathcal{L}^2}$'
    '$||w_\mathrm{A}||_{\mathcal{L}^\infty}$'
    '$||w_\mathrm{B}||_{\mathcal{L}^\infty}$'};
for k=1:4
    figure;
    plot_atlas_2d(runid, 1, 2, k+7);
    xlim(omegaRange); ylim(epsRange);
    xlabel('$\Omega$','FontSize',16,'Interpreter',"latex");
    ylabel('$\epsilon$','FontSize',16,'Interpreter',"latex");
    zlabel(titles{k},'FontSize',16,'Interpreter',"latex");
end
% FRS in reduced coordinates
figure;
plot_atlas_2d(runid, 1, 2, 12);
xlim(omegaRange); ylim(epsRange);
xlabel('$\Omega$','FontSize',16,'Interpreter',"latex");
ylabel('$\epsilon$','FontSize',16,'Interpreter',"latex");
zlabel('$\rho_1$','FontSize',16,'Interpreter',"latex");
figure;
plot_atlas_2d(runid, 1, 2, 13);
xlim(omegaRange); ylim(epsRange);
xlabel('$\Omega$','FontSize',16,'Interpreter',"latex");
ylabel('$\epsilon$','FontSize',16,'Interpreter',"latex");
zlabel('$\rho_2$','FontSize',16,'Interpreter',"latex");